Compute Oxygen Emissions

Part 1: Data Acquisition & Preparation

What we do: 1. Load a hand‐coded World Bank Country→Region mapping. 2. Download per-capita electricity generation mixes from OWID. 3. Download grid carbon intensity (g CO₂e /kWh) from OWID. 4. Merge everything and keep only the most recent year with a non-missing CI per country.

Why: - We need country-level grid carbon intensities and generation shares to scale our oxygen LCA. - Regions are used later for delivery‐mix weighting.

# 1.1 Read region labels (World Bank country & lending groups)
lbl_data <- read_xlsx(here("lbl_regions.xlsx"))
# ‣ Columns: Entity (OWID name), Code (ISO-3), Region (WB group)

# 1.2 Download per-capita electricity mix
grid_data <- read_csv(
  "https://ourworldindata.org/grapher/per-capita-electricity-fossil-nuclear-renewables.csv?v=1&csvType=full&useColumnShortNames=true"
) |>
  rename(
    fossil_pc = per_capita_fossil_generation__kwh_chart_per_capita_electricity_fossil_nuclear_renewables,
    nuclear_pc = per_capita_nuclear_generation__kwh_chart_per_capita_electricity_fossil_nuclear_renewables,
    renewable_pc = per_capita_renewable_generation__kwh_chart_per_capita_electricity_fossil_nuclear_renewables
  )

# 1.3 Download carbon intensity (g CO₂e/kWh)
carbon_intensity <- read_csv(
  "https://ourworldindata.org/grapher/carbon-intensity-electricity.csv?v=1&csvType=full&useColumnShortNames=true"
)

# 1.4 Merge, fallback, filter & compute generation shares
grid_data_ci <- grid_data |>
  # join in Region; if missing, use OWID name as region
  left_join(lbl_data |> select(-Entity), by = "Code") |>
  mutate(Region = if_else(is.na(Region), Entity, Region)) |>
  # attach CO₂ intensity and drop missing
  left_join(
    carbon_intensity |>
      select(Entity, Year, co2_intensity = co2_intensity__gco2_kwh),
    by = c("Entity", "Year")
  ) |>
  drop_na(co2_intensity) |>
  # pick most recent per Entity
  group_by(Entity) |>
  arrange(desc(Year), .by_group = TRUE) |>
  slice(1) |>
  ungroup() |>
  # compute percent shares of generation mix
  mutate(
    total_generation = fossil_pc + nuclear_pc + renewable_pc,
    fossil_generation = 100 * fossil_pc / total_generation,
    nuclear_generation = 100 * nuclear_pc / total_generation,
    renewable_generation = 100 * renewable_pc / total_generation
  ) |>
  # re-attach co2_intensity (to preserve Year) and fix Europe code
  select(-co2_intensity) |>
  left_join(
    carbon_intensity |>
      select(Entity, Year, co2_intensity = co2_intensity__gco2_kwh),
    by = c("Entity", "Year")
  ) |>
  mutate(Code = if_else(Entity == "Europe", "EU", Code))

Part 2: Define Emissions for Each Delivery System

Goal: Break down system‐level CO₂e into legs (“prod”, “transp”, “otro”), scale to 100% purity, and convert to g CO₂e per pure liter.

Key logic: - We measure a total GWP per liter at some purity → allocate by % to each leg → “undo” purity → convert to grams.

# Helper fn: compute component‐level GWP at 100% purity
make_system_emissions <- function(
  components, # leg names
  percent_emissions, # numeric vector summing to 100
  total_kg_GWP_1L, # observed kg CO2e / L at measured purity
  purity_factor # fraction (e.g. 0.995)
) {
  tibble(
    component = components,
    percent_emissions = percent_emissions
  ) |>
    mutate(
      # raw kg per component at measured purity
      kg_GWP_1L_raw = percent_emissions / 100 * total_kg_GWP_1L,
      # scale to 100% purity
      kg_GWP_1L = kg_GWP_1L_raw / purity_factor,
      # grams per pure L
      g_GWP_1L = kg_GWP_1L * 1000
    ) |>
    select(component, percent_emissions, kg_GWP_1L, g_GWP_1L)
}

# System 1: Bulk LOX
system1_lox <- make_system_emissions(
  components = c("total", "prod", "transp", "otro"),
  percent_emissions = c(100, 82, 18, 0),
  total_kg_GWP_1L = 0.0001700,
  purity_factor = 0.995
)

# System 2: Cylinder LOX
system2_cylox <- make_system_emissions(
  components = c("total", "prod", "transp", "otro"),
  percent_emissions = c(100, 34, 36, 30),
  total_kg_GWP_1L = 0.0004110,
  purity_factor = 0.995
)

# System 3: PSA_local
system3_psal <- make_system_emissions(
  components = c("total", "prod", "transp", "otro"),
  percent_emissions = c(100, 99, 0, 1),
  total_kg_GWP_1L = 0.0000844,
  purity_factor = 0.93
)

# System 4: Oxygen concentrator
system4_oxconc <- make_system_emissions(
  components = c("total", "prod", "transp", "otro"),
  percent_emissions = c(100, 99.5, 0, 0.5),
  total_kg_GWP_1L = 0.000123,
  purity_factor = 0.93
)

# System 5: Composite PSA cylinder
#  - prod & otro from PSA local (93%)
#  - transp from cylinder LOX (undo 99.5%→93%)
purity_cylox <- 0.995
purity_psal <- 0.93

psal_base <- system3_psal |>
  filter(component %in% c("prod", "otro")) |>
  select(component, kg_GWP_1L, g_GWP_1L)

cylox_transp_adj <- system2_cylox |>
  filter(component == "transp") |>
  mutate(
    raw_kg = kg_GWP_1L * purity_cylox,
    kg_GWP_1L = raw_kg / purity_psal,
    g_GWP_1L = kg_GWP_1L * 1000
  ) |>
  select(component, kg_GWP_1L, g_GWP_1L)

components_only <- bind_rows(psal_base, cylox_transp_adj)
new_total_kg <- sum(components_only$kg_GWP_1L)
new_total_g <- sum(components_only$g_GWP_1L)

system5_psac <- bind_rows(
  tibble(
    component = "total",
    percent_emissions = 100,
    kg_GWP_1L = new_total_kg,
    g_GWP_1L = new_total_g
  ),
  components_only |>
    mutate(percent_emissions = 100 * kg_GWP_1L / new_total_kg)
) |>
  arrange(factor(component, levels = c("total", "prod", "transp", "otro")))

Liquid Oxygen

Cylinder of LOX

PSA Piped

Oxygen Concentrator

PSA Cylinder

Part 3: Calculate Emissions per Country per Delivery System

Objective: For each country and each of our five delivery systems, compute the total CO₂e (g) per pure-L of O₂ by: 1. Combining the per-system component tables into one long format. 2. Cross-joining with each country’s grid carbon intensity (CI) and scaling the “production” leg proportionally to the country’s CI. 3. Re-normalizing the component shares to sum to 100%, summing the grams CO₂e across legs, pivoting to a wide layout, and attaching these values back to the country CI table.

# 3.1 Combine all five system breakdowns into one long table
#     - 'systems' is a named list of tibbles: one per delivery system
#     - imap_dfr() binds them row-wise and adds a 'system' column
systems <- list(
  system1_lox = system1_lox,
  system2_cylox = system2_cylox,
  system3_psal = system3_psal,
  system4_oxconc = system4_oxconc,
  system5_psac = system5_psac
)

systems_long <- imap_dfr(
  systems,
  ~ mutate(.x, system = .y) # .x = each tibble; .y = its name
)

# 3.2 Cross-join each country with every system component,
#     then scale only the “prod” (production) leg according to the ratio:
#       country_CI / ONTARIO_BASELINE_CI
ONTARIO_BASELINE_CI <- 28 # g CO₂e per kWh from Ontario LCA

emissions_components <- grid_data_ci %>%
  select(Entity, co2_intensity) %>% # keep only country name & its CI
  crossing(systems_long) %>% # replicate each component for each country
  mutate(
    # if this row is the production leg, scale g_GWP_1L by CI ratio
    g_GWP_1L = if_else(
      component == "prod",
      g_GWP_1L * (co2_intensity / ONTARIO_BASELINE_CI),
      g_GWP_1L
    )
  ) %>%
  group_by(Entity, system) %>%
  mutate(
    # re-compute each leg’s % share so that prod + transp + otro = 100
    percent_emissions = 100 * g_GWP_1L / sum(g_GWP_1L[component != "total"])
  ) %>%
  ungroup()

# 3.3 Summarize to one total g CO₂e per pure-L for each country × system,
#     then pivot to wide format and join back to the country CI table.
emissions_totals <- emissions_components %>%
  filter(component != "total") %>% # drop the original 'total' row
  group_by(Entity, system) %>%
  summarise(
    total_g_GWP_1L = sum(g_GWP_1L), # grams CO₂e per pure-L
    .groups = "drop"
  )

emissions_wide <- emissions_totals %>%
  # prefix system names with "GWP_" for clarity in column names
  mutate(system = paste0("GWP_", system)) %>%
  pivot_wider(
    names_from = system,
    values_from = total_g_GWP_1L
  )

# Attach the per-system GWP columns to the original country CI table
country_emissions <- grid_data_ci %>%
  left_join(emissions_wide, by = "Entity")

Part 4: Calculate emission 1 L/min per country

Objective Blend the five delivery‐system GWP factors with region‐specific usage shares to derive a single g CO₂e per pure-L of O₂ per minute for every country.

Key Assumptions

• Regional Shares The oxygen_distribution table reflects typical delivery-system splits by region (Supplemental Table X). Shares sum to 100% in each region.
• Continental Aggregation Rows labeled “Africa”, “Europe”, etc., are not raw data but aggregated averages of their constituent sub-regions.

Steps

1. Define regional delivery-system mix shares (Supplemental Table X).
2. Compute weighted average g CO₂e per pure-L using those shares.
3. Compute continent-level means for “whole” regions and override.
4. Save oxemissions_bycountry.csv.

4.1 Delivery-mix tribble

# 4.2 Prepare the regional mix shares for weighting ---------------------------

# Append “_share” to each system column so it lines up with the GWP columns
oxygen_mix <- oxygen_distribution %>%
  rename_with(~ paste0(.x, "_share"), -Entity)

# List of the per-system GWP columns produced in Part 3
system_cols <- c(
  "GWP_system1_lox",
  "GWP_system2_cylox",
  "GWP_system3_psal",
  "GWP_system4_oxconc",
  "GWP_system5_psac"
)

# Corresponding list of the share columns from oxygen_mix
share_cols <- c(
  "system1_lox_share",
  "system2_cylox_share",
  "system3_psal_share",
  "system4_oxconc_share",
  "system5_psac_share"
)


# 4.3 Compute weighted average emissions per pure-L for each country -----------

country_1l <- country_emissions %>%
  # bring in the regional system shares
  left_join(oxygen_mix, by = c("Region" = "Entity")) %>%
  # rowwise allows us to compute a per-row sum across multiple columns
  rowwise() %>%
  mutate(
    # sum(GWP_i * share_i) for i = 1…5, then divide by 100 because shares are %
    g_CO2e_1L = sum(
      c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
      na.rm = TRUE
    ) /
      100
  ) %>%
  ungroup()


# 4.4 Compute continent-level averages and override continent rows ------------

# Which Entity names represent whole continents
continent_names <- c("Africa", "Asia", "Europe", "Oceania", "South America")

# 1) Exclude continent rows, map each sub-region to its continent,
#    then compute the mean g_CO2e_1L for each continent
continent_avgs <- country_1l %>%
  filter(!Entity %in% continent_names) %>%
  mutate(
    continent = case_when(
      Region %in% c("SSA", "MENA") ~ "Africa",
      str_detect(Region, regex("Asia", ignore_case = TRUE)) ~ "Asia",
      str_detect(Region, regex("Europe", ignore_case = TRUE)) ~ "Europe",
      str_detect(Region, regex("Oceania", ignore_case = TRUE)) ~ "Oceania",
      str_detect(Region, regex("Latin America", ignore_case = TRUE)) ~
        "South America",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(continent)) %>%
  group_by(continent) %>%
  summarise(
    g_CO2e_1L = mean(g_CO2e_1L, na.rm = TRUE),
    .groups = "drop"
  )

# 2) Join those continent means back to the full table
#    and overwrite the continent rows with their average
country_final <- country_1l %>%
  left_join(continent_avgs, by = c("Entity" = "continent")) %>%
  mutate(
    g_CO2e_1L = if_else(
      Entity %in% continent_names,
      g_CO2e_1L.y, # continent mean
      g_CO2e_1L.x # original country value
    )
  ) %>%
  # drop the helper columns and any zero‐emission rows
  select(-g_CO2e_1L.x, -g_CO2e_1L.y) %>%
  filter(g_CO2e_1L > 0)


# 4.5 Render the final table as an interactive DataTable -----------------------

DT::datatable(
  country_final,
  extensions = 'Buttons',
  options = list(
    pageLength = 25, # show 25 rows at a time
    scrollX = TRUE, # horizontal scrolling
    dom = 'Bfrtip', # show Buttons, filter, table, info, pagination
    buttons = c('copy', 'csv', 'excel', 'pdf')
  ),
  caption = '### Baseline O₂ Emissions per L by Country'
)

Part 5: Projected Grid Carbon Intensity Scenarios for 2030

In this section we augment our baseline country data with three hypothetical future grid carbon‐intensity metrics (g CO₂e /kWh), under differing decarbonization pathways:

1. 2030 Renewable Target Scenario (co2_intensity_2030)
   • Assumes each country achieves its stated 2030 renewable energy share
   • Nuclear share remains fixed at today’s level
   • Fossil reduction is entirely backfilled by renewables
2. Half-Fossil Scenario (co2_intensity_half)
   • Assumes current fossil generation share is cut in half
   • Nuclear share remains constant
   • Displaced fossil generation replaced by renewables
3. Clean-Grid Scenario (co2_intensity_clean)
   • Assumes 100 % elimination of fossil generation
   • Nuclear share remains at baseline
   • Residual energy supplied solely by renewables

---

Key Assumptions & Parameters

• Nuclear emission factor: 12 g CO₂e /kWh (IPCC AR5)
• Pure-fossil & pure-renewable factors: Median carbon intensities of grids with > 99 % fossil or > 99 % renewables
• Imputation Rules for Missing Targets:
  • European countries with no target → use “Europe” (EU) average
  • South America → inherits mean of all Latin America sub-regions
  • Remaining Latin America gaps → filled with the South America value

We compute these three new intensity columns in our target_intensities table, then feed them into Part 6’s per-liter emissions calculations.```

# 5.1 Read & merge 2030 RES targets
targets_renwbl <- read_xlsx(here("targets_download.xlsx"), sheet=4) |>
  select(Code=country_code, target_perc_2030=res_share_target) |>
  mutate(Code=if_else(Code=="XKX","OWID_KOS",Code))

target_2030 <- country_final |>
  left_join(targets_renwbl, by="Code") |>
  mutate(
    target_perc_2030 = if_else(
      is.na(target_perc_2030) & str_starts(Region,"Europe"),
      first(target_perc_2030[Entity=="Europe"], NA_real_),
      target_perc_2030
    ),
    target_perc_2030 = if_else(
      Entity=="South America",
      mean(target_perc_2030[str_starts(Region,"Latin America")], na.rm=TRUE),
      target_perc_2030
    ),
    target_perc_2030 = if_else(
      is.na(target_perc_2030) & str_starts(Region,"Latin America"),
      first(target_perc_2030[Entity=="South America"], NA_real_),
      target_perc_2030
    )
  ) |>
  filter(!is.na(Code))

# 5.2 Define pure emission factors
ef_nuclear <- 12
ef_fossil  <- grid_data_ci |> filter(fossil_generation>99)   |> summarise(median(co2_intensity,na.rm=TRUE)) |> pull()
ef_renew   <- grid_data_ci |> filter(renewable_generation>99)|> summarise(median(co2_intensity,na.rm=TRUE)) |> pull()

# 5.3 Compute scenario intensities
target_intensities <- target_2030 |>
  mutate(
    co2_intensity_2030 = if_else(
      is.na(target_perc_2030), NA_real_,
      {
        rg <- pmin(100 - nuclear_generation, pmax(renewable_generation, target_perc_2030))
        fg <- 100 - nuclear_generation - rg
        fg/100*ef_fossil + rg/100*ef_renew + nuclear_generation/100*ef_nuclear
      }
    ),
    co2_intensity_half = {
      fg <- fossil_generation * 0.5
      rg <- 100 - nuclear_generation - fg
      fg/100*ef_fossil + rg/100*ef_renew + nuclear_generation/100*ef_nuclear
    },
    co2_intensity_clean = {
      rg <- 100 - nuclear_generation
      rg/100*ef_renew + nuclear_generation/100*ef_nuclear
    }
  )

Part 6: Apply Scenario Intensities & Compile Final Emissions Table

In this final stage, we take the three grid-intensity projections from Part 5 and feed them through our delivery-system LCA pipeline (as in Part 3), apply the regional mixes from Part 4, and then stitch everything together into one comprehensive dataset.

Workflow Overview: 1. Inject Scenario CI - For each scenario (2030 target, half-fossil, clean-grid), replace the baseline co2_intensity in grid_data_ci with the scenario-specific column (co2_intensity_2030, co2_intensity_half, or co2_intensity_clean). 2. Recompute System-Level Emissions - Cross-join with systems_long to get every country × delivery-system component. - Scale the production leg by (scenario CI / Ontario baseline CI). - Recalculate component shares so that prod + transp + other = 100 %. 3. Aggregate to Country × System - Sum the grams CO₂e per pure-L for each system, then pivot to wide format (e.g. GWP_system1_lox, …). 4. Compute 1 L/min Emissions per Country - Apply the regional delivery-mix shares to each system’s wide table to get a weighted g_CO2e_1L_<scenario> for every country. 5. Merge Baseline & Projections - Left-join baseline and all three scenario-specific g_CO2e_1L columns into a single final_table. - Render as an interactive datatable or export to CSV.

Final Output (final_table): - Identifiers: Entity, Code, Region - Baseline Metrics: co2_intensity, g_CO2e_1L - 2030 Target Scenario: co2_intensity_2030, g_CO2e_1L_2030 - Half-Fossil Scenario: co2_intensity_half, g_CO2e_1L_half - Clean-Grid Scenario: co2_intensity_clean, g_CO2e_1L_clean

# 6.1 2030 target
grid_2030 <- grid_data_ci |>
  select(-co2_intensity) |>
  distinct(Code, .keep_all = TRUE) |>
  left_join(
    select(target_intensities, Code, co2_intensity_2030),
    by = "Code"
  ) |>
  rename(co2_intensity = co2_intensity_2030)

emissions_components_2030 <- grid_2030 |>
  select(Entity = Code, co2_intensity) |>
  crossing(systems_long) |>
  mutate(
    g_GWP_1L = if_else(
      component == "prod",
      g_GWP_1L * (co2_intensity / 28),
      g_GWP_1L
    )
  ) |>
  group_by(Entity, system) |>
  mutate(
    percent_emissions = g_GWP_1L / sum(g_GWP_1L[component != "total"]) * 100
  ) |>
  ungroup()

emissions_totals_2030 <- emissions_components_2030 |>
  filter(component != "total") |>
  group_by(Entity, system) |>
  summarise(total_g_GWP_1L = sum(g_GWP_1L), .groups = "drop")

emissions_wide_2030 <- emissions_totals_2030 |>
  mutate(system = paste0("GWP_", system)) |>
  pivot_wider(names_from = system, values_from = total_g_GWP_1L)

country_emissions_2030 <- grid_2030 |>
  select(Code) |>
  left_join(emissions_wide_2030, by = c("Code" = "Entity")) |>
  filter(!is.na(Code) & Code != "")

country_1l_2030 <- country_emissions_2030 |>
  left_join(select(grid_2030, Code, Region), by = "Code") |>
  left_join(oxygen_mix, by = c("Region" = "Entity")) |>
  rowwise() |>
  mutate(
    all_gwp_na = if_all(all_of(system_cols), is.na),
    g_CO2e_1L_2030 = if_else(
      all_gwp_na,
      NA_real_,
      sum(
        c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
        na.rm = TRUE
      ) /
        100
    )
  ) |>
  ungroup() |>
  select(Code, g_CO2e_1L_2030) |>
  filter(!is.na(g_CO2e_1L_2030))

# 6.2 –50% fossil
grid_half <- grid_data_ci |>
  select(-co2_intensity) |>
  distinct(Code, .keep_all = TRUE) |>
  left_join(
    select(target_intensities, Code, co2_intensity_half),
    by = "Code"
  ) |>
  rename(co2_intensity = co2_intensity_half)

emissions_components_half <- grid_half |>
  select(Entity = Code, co2_intensity) |>
  crossing(systems_long) |>
  mutate(
    g_GWP_1L = if_else(
      component == "prod",
      g_GWP_1L * (co2_intensity / 28),
      g_GWP_1L
    )
  ) |>
  group_by(Entity, system) |>
  mutate(
    percent_emissions = g_GWP_1L / sum(g_GWP_1L[component != "total"]) * 100
  ) |>
  ungroup()

emissions_totals_half <- emissions_components_half |>
  filter(component != "total") |>
  group_by(Entity, system) |>
  summarise(total_g_GWP_1L = sum(g_GWP_1L), .groups = "drop")

emissions_wide_half <- emissions_totals_half |>
  mutate(system = paste0("GWP_", system)) |>
  pivot_wider(names_from = system, values_from = total_g_GWP_1L)

country_emissions_half <- grid_half |>
  select(Code) |>
  left_join(emissions_wide_half, by = c("Code" = "Entity")) |>
  filter(!is.na(Code) & Code != "")

country_1l_half <- country_emissions_half |>
  left_join(select(grid_half, Code, Region), by = "Code") |>
  left_join(oxygen_mix, by = c("Region" = "Entity")) |>
  rowwise() |>
  mutate(
    all_gwp_na = if_all(all_of(system_cols), is.na),
    g_CO2e_1L_half = if_else(
      all_gwp_na,
      NA_real_,
      sum(
        c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
        na.rm = TRUE
      ) /
        100
    )
  ) |>
  ungroup() |>
  select(Code, g_CO2e_1L_half) |>
  filter(!is.na(g_CO2e_1L_half))

# 6.3 0% fossil
grid_clean <- grid_data_ci |>
  select(-co2_intensity) |>
  distinct(Code, .keep_all = TRUE) |>
  left_join(
    select(target_intensities, Code, co2_intensity_clean),
    by = "Code"
  ) |>
  rename(co2_intensity = co2_intensity_clean)

emissions_components_clean <- grid_clean |>
  select(Entity = Code, co2_intensity) |>
  crossing(systems_long) |>
  mutate(
    g_GWP_1L = if_else(
      component == "prod",
      g_GWP_1L * (co2_intensity / 28),
      g_GWP_1L
    )
  ) |>
  group_by(Entity, system) |>
  mutate(
    percent_emissions = g_GWP_1L / sum(g_GWP_1L[component != "total"]) * 100
  ) |>
  ungroup()

emissions_totals_clean <- emissions_components_clean |>
  filter(component != "total") |>
  group_by(Entity, system) |>
  summarise(total_g_GWP_1L = sum(g_GWP_1L), .groups = "drop")

emissions_wide_clean <- emissions_totals_clean |>
  mutate(system = paste0("GWP_", system)) |>
  pivot_wider(names_from = system, values_from = total_g_GWP_1L)

country_emissions_clean <- grid_clean |>
  select(Code) |>
  left_join(emissions_wide_clean, by = c("Code" = "Entity")) |>
  filter(!is.na(Code) & Code != "")

country_1l_clean <- country_emissions_clean |>
  left_join(select(grid_clean, Code, Region), by = "Code") |>
  left_join(oxygen_mix, by = c("Region" = "Entity")) |>
  rowwise() |>
  mutate(
    all_gwp_na = if_all(all_of(system_cols), is.na),
    g_CO2e_1L_clean = if_else(
      all_gwp_na,
      NA_real_,
      sum(
        c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
        na.rm = TRUE
      ) /
        100
    )
  ) |>
  ungroup() |>
  select(Code, g_CO2e_1L_clean) |>
  filter(!is.na(g_CO2e_1L_clean))

# 6.4 Final assembly & export
final_table <- target_intensities |>
  select(
    Entity,
    Code,
    Region,
    co2_intensity,
    co2_intensity_2030,
    co2_intensity_half,
    co2_intensity_clean,
    g_CO2e_1L
  ) |>
  left_join(country_1l_2030, by = "Code") |>
  left_join(country_1l_half, by = "Code") |>
  left_join(country_1l_clean, by = "Code")


# ── 6.5 Render final_table as an interactive DataTable ────────────────────────

DT::datatable(
  final_table,
  extensions = 'Buttons',
  options = list(
    pageLength = 25,
    scrollX = TRUE,
    dom = 'Bfrtip',
    buttons = c('csv', 'excel')
  ),
  caption = '### Emissions per Liter Projections'
)